//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   DD_1D_Equations.i.hh
 * \author Jeremy Roberts
 * \date   Sep 4, 2011
 * \brief  DD_1D_Equations inline member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef DD_1D_EQUATIONS_I_HH_
#define DD_1D_EQUATIONS_I_HH_

namespace slabtran
{

// Solve in a cell.
template <class MF, class PCF>
inline void DD_1D_Equations::solve(const int                  i,
                                   const double               sigma,
                                   const double               q,
                                   const Omega                &mu,
                                   const double               wt,
                                   const Moment_to_Discrete_t &M,
                                   const Space_Vector         &in_psi,
                                   MF                         &phi,
                                   Space_Vector               &out_psi,
                                   PCF                        &psi,
                                   const int                  angle,
                                   bool                       store_psi)
{
  // The absolute value on the angle lets us use the same
  // code for both directions.
  double d_factor = 2.0 * fabs(mu) * d_mesh->inv_width(i, I);

  // Calculate cell-centered angular flux.
  d_psi = (q + d_factor * in_psi) / (sigma + d_factor);

  // If requested, store the cell-centered angular flux.
  if (store_psi) psi(i, 0, angle) = d_psi;

  // Calculate new exiting flux.  This is the DD definition!
  out_psi = 2.0 * d_psi - in_psi;

  // Loop through flux moments and add up the contributions.  We make use
  // of the cardinal index interfaces, as i and l are the cardinal indices
  // of cells and moments in 1-d.
  for (int l = 0, N = phi.legendre_order(); l <= N; l++)
    phi(i, 0, l) += d_psi * wt * M(angle, l);
  /// \todo [PERFORMANCE] Use BLAS call here for loop (need direct access to
  ///       continuous data)
}

} // end namespace slabtran

#endif /* DD_1D_EQUATIONS_I_HH_ */

//---------------------------------------------------------------------------//
//              end of DD_1D_Equations.i.hh
//---------------------------------------------------------------------------//
